/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "PatchToPatchInterpolationTemplate.H"
#include "objectHit.H"
#include "pointHit.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class FromPatch, class ToPatch>
void PatchToPatchInterpolation<FromPatch, ToPatch>::calcPointAddressing() const
{
	// Calculate pointWeights

	pointWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.nPoints());
	FieldField<Field, scalar>& pointWeights = *pointWeightsPtr_;

	pointDistancePtr_ = new scalarField(toPatch_.nPoints(), GREAT);
	scalarField& pointDistance = *pointDistancePtr_;

	const pointField& fromPatchPoints = fromPatch_.localPoints();
	const List<typename FromPatch::FaceType>& fromPatchFaces =
		fromPatch_.localFaces();

	const pointField& toPatchPoints = toPatch_.localPoints();
	const vectorField& projectionDirection = toPatch_.pointNormals();
	const edgeList& toPatchEdges = toPatch_.edges();
	const labelListList& toPatchPointEdges = toPatch_.pointEdges();

	if (debug)
	{
		Info << "projecting points" << endl;
	}

	List<objectHit> proj =
		toPatch_.projectPoints(fromPatch_, projectionDirection, alg_, dir_);

	pointAddressingPtr_ = new labelList(proj.size(), -1);
	labelList& pointAddressing = *pointAddressingPtr_;

	bool doWeights = false;

	forAll (pointAddressing, pointI)
	{
		doWeights = false;

		const typename FromPatch::FaceType& hitFace =
			fromPatchFaces[proj[pointI].hitObject()];

		point hitPoint = point::zero;

		if (proj[pointI].hit())
		{
			// A hit exists
			doWeights = true;

			pointAddressing[pointI] = proj[pointI].hitObject();

			pointHit curHit =
				hitFace.ray
				(
					toPatchPoints[pointI],
					projectionDirection[pointI],
					fromPatchPoints,
					alg_,
					dir_
				);

			// Grab distance to target
			if (dir_ == intersection::CONTACT_SPHERE)
			{
				pointDistance[pointI] =
					hitFace.contactSphereDiameter
					(
					    toPatchPoints[pointI],
					    projectionDirection[pointI],
					    fromPatchPoints
					);
			}
			else
			{
				pointDistance[pointI] = curHit.distance();
			}

			// Grab hit point
			hitPoint = curHit.hitPoint();
		}
		else if (projectionTol_() > SMALL)
		{
			// Check for a near miss
			pointHit ph =
				hitFace.ray
				(
					toPatchPoints[pointI],
					projectionDirection[pointI],
					fromPatchPoints,
					alg_,
					dir_
				);

			scalar dist =
				Foam::mag
				(
					toPatchPoints[pointI]
				  + projectionDirection[pointI]*ph.distance()
				  - ph.missPoint()
				);

			// Calculate the local tolerance
			scalar minEdgeLength = GREAT;

			// Do shortest edge of hit object
			edgeList hitFaceEdges =
				fromPatchFaces[proj[pointI].hitObject()].edges();

			forAll (hitFaceEdges, edgeI)
			{
				minEdgeLength =
					min
					(
					    minEdgeLength,
					    hitFaceEdges[edgeI].mag(fromPatchPoints)
					);
			}

			const labelList& curEdges = toPatchPointEdges[pointI];

			forAll (curEdges, edgeI)
			{
				minEdgeLength =
					min
					(
					    minEdgeLength,
					    toPatchEdges[curEdges[edgeI]].mag(toPatchPoints)
					);
			}

			if (dist < minEdgeLength*projectionTol_())
			{
				// This point is being corrected
				doWeights = true;

				pointAddressing[pointI] = proj[pointI].hitObject();

				// Grab nearest point on face as hit point
				hitPoint = ph.missPoint();

				// Grab distance to target
				if (dir_ == intersection::CONTACT_SPHERE)
				{
					pointDistance[pointI] =
					    hitFace.contactSphereDiameter
					    (
					        toPatchPoints[pointI],
					        projectionDirection[pointI],
					        fromPatchPoints
					    );
				}
				else
				{
					pointDistance[pointI] =
					    (
					        projectionDirection[pointI]
					        /mag(projectionDirection[pointI])
					    )
					  & (hitPoint - toPatchPoints[pointI]);
				}
			}
		}

		if (doWeights)
		{
			// Set interpolation pointWeights
			pointWeights.set(pointI, new scalarField(hitFace.size()));

			pointField hitFacePoints = hitFace.points(fromPatchPoints);

			forAll (hitFacePoints, masterPointI)
			{
				pointWeights[pointI][masterPointI] =
					1.0/
					(
					    mag
					    (
					        hitFacePoints[masterPointI]
					      - hitPoint
					    )
					  + VSMALL
					);
			}

			pointWeights[pointI] /= sum(pointWeights[pointI]);
		}
		else
		{
			pointWeights.set(pointI, new scalarField(0));
		}
	}
}


template<class FromPatch, class ToPatch>
void PatchToPatchInterpolation<FromPatch, ToPatch>::calcFaceAddressing() const
{
	faceWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.size());
	FieldField<Field, scalar>& faceWeights = *faceWeightsPtr_;

	faceDistancePtr_ = new scalarField(toPatch_.size(), GREAT);
	scalarField& faceDistance = *faceDistancePtr_;

	if (debug)
	{
		Info << "projecting face centres" << endl;
	}

	const pointField& fromPatchPoints = fromPatch_.points();
	const typename FromPatch::FaceListType& fromPatchFaces = fromPatch_;
	const labelListList& fromPatchFaceFaces = fromPatch_.faceFaces();

	vectorField fromPatchFaceCentres(fromPatchFaces.size());

	forAll (fromPatchFaceCentres, faceI)
	{
		fromPatchFaceCentres[faceI] =
			fromPatchFaces[faceI].centre(fromPatchPoints);
	}

	const pointField& toPatchPoints = toPatch_.points();
	const typename ToPatch::FaceListType& toPatchFaces = toPatch_;

	const vectorField& projectionDirection = toPatch_.faceNormals();

	List<objectHit> proj =
		toPatch_.projectFaceCentres
		(
			fromPatch_,
			projectionDirection,
			alg_,
			dir_
		);

	faceAddressingPtr_ = new labelList(proj.size(), -1);
	labelList& faceAddressing = *faceAddressingPtr_;

	forAll (faceAddressing, faceI)
	{
		if (proj[faceI].hit())
		{
			// A hit exists
			faceAddressing[faceI] = proj[faceI].hitObject();

			const typename FromPatch::FaceType& hitFace =
				fromPatchFaces[faceAddressing[faceI]];

			pointHit curHit =
				hitFace.ray
				(
					toPatchFaces[faceI].centre(toPatchPoints),
					projectionDirection[faceI],
					fromPatchPoints,
					alg_,
					dir_
				);

			// grab distance to target
			faceDistance[faceI] = curHit.distance();

			// grab face centre of the hit face
			const point& hitFaceCentre =
				fromPatchFaceCentres[faceAddressing[faceI]];

			// grab neighbours of hit face
			const labelList& neighbours =
				fromPatchFaceFaces[faceAddressing[faceI]];

			scalar m = mag(curHit.hitPoint() - hitFaceCentre);

			if
			(
				m < directHitTol_()                           // Direct hit
			 || neighbours.empty()
			)
			{
				faceWeights.set(faceI, new scalarField(1));
				faceWeights[faceI][0] = 1.0;
			}
			else
			{
				// set interpolation faceWeights

				// The first coefficient corresponds to the centre face.
				// The rest is ordered in the same way as the faceFaces list.
				faceWeights.set(faceI, new scalarField(neighbours.size() + 1));

				faceWeights[faceI][0] = 1.0/m;

				forAll (neighbours, nI)
				{
					faceWeights[faceI][nI + 1] =
					1.0/
					(
					    mag
					    (
					        fromPatchFaceCentres[neighbours[nI]]
					      - curHit.hitPoint()
					    )
					  + VSMALL
					);
				}
			}

			faceWeights[faceI] /= sum(faceWeights[faceI]);
		}
		else
		{
			faceWeights.set(faceI, new scalarField(0));
		}
	}
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //
